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We present a combination of tools which allows for investigation of the coupled orbital 
and rotational dynamics of two rigid bodies with nearly arbitrary shape and mass distribu- 
tion, under the influence of their mutual gravitational potential. Methods for calculating 
that mutual potential and resulting forces and moments for a polyhedral body represen- 
tation are simple and efficient. Discrete equations of motion, referred to as the Lie Group 
Variational Integrator (LGVI), preserve the structure of the configuration space, SE(3), as 
well as the geometric features represented by the total energy and the total angular mo- 
mentum. The synthesis of these approaches allows us to simulate the full two body problem 
accurately and efficiently. Simulation results are given for two octahedral rigid bodies for 
comparison with other integration methods and to show the qualities of the results thus 
obtained. A significant improvement is seen over other integration methods while correctly 
capturing the interesting effects of strong orbit and attitude dynamics coupling, in multiple 
scenarios. 

I. Introduction 

The full two body problem studies the dynamics of two irregular rigid bodies interacting under a mutual 
potential. The full two body problem arises in numerous engineering and scientific fields. Our focus is on 
the dynamics of two rigid bodies in space due to their mutual gravitational potential. This depends on both 
the relative position and the relative attitude of the bodies. Therefore, the translational orbit dynamics and 
the rotational attitude dynamics are coupled in the full two body problem. For example, the trajectory of a 
very large spacecraft around the Earth is affected by the attitude of the spacecraft, and the dynamics of a 
binary asteroid pair are characterized by the non-spherical mass distributions of the two bodies. Recently, 
interest in the full (two) body problem as applied to binary objects in space has increased, due to evidence 
that such binary systems are common among the overall asteroid population and form an especially large 
percentage of near-earth objects, with some estimates as high as 16 percent. 1-4 In addition, advances in 
radar and optical observation methods 5 have allowed for better modelling of the shapes and characteristics 
of bodies in binary systems, encouraging further study of their dynamics. 

Some previous work of relevance includes Maciejewski's presentation of equations of motion of the full 
two body problem in incrtial and relative coordinates. 6 He also discussed the existence of relative equilibria. 
Scheeres derived a stability condition for the full two body problem, 7 and he studied the planar stability of 
an ellipsoid-sphere model. 8 Spacecraft motion about binary asteroids has been discussed using the restricted 
three body model, 9 ' 10 and the four body model. 11 

The mutual gravitational potential of two rigid celestial bodies has been expressed using spherical har- 
monics. 12 ' 13 But, the harmonic expansion is not guaranteed to converge. Convergence is shown to be an 
unstable property of such spherical harmonic series. 14 By this we mean that an arbitrarily small change to 
the mass distribution may cause a previously convergent series to diverge. Another commonly used approach 
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for evaluating the mutual gravitational potential is to fill each rigid body's volume with a distribution of point 
masses, fixed with respect to one another, the sum of which equals the respective body's total mass. 15,16 
Although the mutual potential obtained for two rigid bodies using this approach converges to the true gravity 
field in the limit as the number of point masses becomes arbitrarily large, there are significant errors in the 
computation of gravitational forces from that mutual potential. 17 

These problems are avoided with the approach presented in Ref. 18 for calculating mutual gravitational 
potential based on polyhedral body models. This method converges at all points exterior to the two bodies. 
By representing the bodies as polyhedra, the flexibility of the formalism over more specialized representations 
such as spheres and ellipsoids is increased. A polyhedral rigid body model can directly include important 
features such as craters, caves, or deep clefts where contact binaries meet. The entire body does not have to 
be modelled uniformly at a high resolution. The errors in the mutual potential computation can be reduced 
to the level of error in each body's shape determination, and to the level of discretization chosen for that 
shape. Derivatives of this polyhedral mutual potential formulation are given in Ref. 19 and, with some 
corrections and improvement of notation, in Ref. 20. These derivatives determine the forces and torques 
exerted by the bodies on each other, for use in cither inertial or relative equations of motion that describe 
the full body dynamics. 

These full body dynamics arise from Lagrangian and Hamiltonian mechanics; they are characterized by 
symplectic, momentum and energy preserving properties. These geometric features determine the qualitative 
behavior of the full body dynamics, and they can serve as a basis for further theoretical study of the full 
body problem. The configuration space of the full body dynamics have a Lie group structure referred to 
as the Euclidean group, SE(3). However, general numerical integration methods, including the widely used 
Runge-Kutta schemes, neither preserve the Lie group structure nor these geometric properties. 21 

The variational approach 22 and Lie group methods 23 provide systematic methods of constructing struc- 
ture preserving numerical integrators. The idea of the variational approach is to discretize Hamilton's 
principle rather than the continuous equations of motion. 22 The numerical integrator obtained from the 
discrete Hamilton's principle exhibits excellent energy properties, conserves first integrals, and preserves 
the symplectic structure. Lie group methods consist of numerical integrators that preserve the geometry of 
the configuration space by automatically remaining on the Lie group. 23 A Lie group method is explicitly 
adopted for the variational integrator in Ref. 24 and 25. This unified integrator, hereafter referred to as 
the Lie Group Variational Integrator (or LGVI for short), is symplectic and momentum preserving, and it 
exhibits good total energy behavior for exponentially long time periods. It also preserves the Euclidian Lie 
group structure without the use of local charts, reprojection, or constraints. 

Numerical simulation of the full body problem involves two major problems; a large computational 
burden in computing mutual gravitational forces and moments, and inaccuracy of numerical integrators. 
The forces and moments must constantly be reevaluated with any position change or any orientation change, 
not only at each time step but at each sub-step involved in the differencing scheme behind any general 
numerical integrator. Therefore such general numerical integrators amplify the computation cost for finding 
the forces and moments. The accuracy of such general integrators also rapidly degrades as the simulation 
time increases. They fail to preserve the conserved quantities such as total energy and angular momentum, 
which determine the qualitative behavior of the full body dynamics. Attitude errors tend to accumulate as a 
consequence of numerical errors, and this attitude degradation causes significant errors in the gravitational 
force and moment computation. 

The unified treatment given in this paper of the polyhedral mutual potential formulation combined with 
the LGVI presents a solution to these two major dynamic simulation problems. Using polyhedral models, 
we can approximate irregular bodies to a specified accuracy, and we can control the computational burden 
to compute the mutual gravitational forces and moments by choosing the level of body discretization and 
the number of scries terms employed in our formulation. The LGVI, as presented in this paper, preserves 
the conserved quantities of the full body dynamics as well as the orthogonal structure of the rotation 
matrices. The obtained simulation results exhibit good stability properties for the invariants of motion for 
exponentially long times. The computational load is further minimized since the LGVI requires one force 
and torque evaluation per integration step for second order accuracy. 

Subsequent sections of this paper arc organized as follows. Algorithms for the mutual gravitational forces 
and moments computation for polyhedral body models are presented in section |n] Some description of the 
full two body problem and its mathematical properties and the continuous relative equations of motion are 
given in section IIIII The Lie Group Variational Integrator for computing the dynamics of the full body 
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problem is given in section Hvl Selected simulation results are given in section [V] for two octahedral rigid 
bodies in several scenarios, for comparison with other integration methods and to show the qualities of 
the results thus obtained. Conclusions drawn from the results for these scenarios about the accuracy and 
suitability of our methods, and about computational burden, are given in the last section. 

II. Polyhedral Mutual Gravitational Potential, Forces, and Moments 

In this section, we present computational algorithms for determining the mutual gravitational potential, 
forces, and moments given polyhedral models of each of the two rigid bodies. The algorithm to compute the 
mutual potential relatively efficiently with such a modelling approach is outlined in Ref. 18. The methods 
outlined in Ref. 19, and in Ref. 20 with some corrections and improvement of notation, can then be used to 
compute the gradients of the mutual potential (i.e. forces and moments). A more detailed description and 
derivation of materials in this section, including extensions to computing the forces and moments for incrtial 
equations of motion, can be found in Ref. 20. 

A. Context for Methodology 

In Ref. 18, a uniform density or homogenous mass distribution within the entire volume of each rigid body is 
assumed. This may be a realistic assumption for binary asteroids based on some empirical data. 26 However, 
the underlying method does also allow for approximating arbitrary density variations, as required for greater 
asteroid modelling fidelity and for systems in which one or more bodies is a spacecraft. This is mentioned 
without proper explanation in Ref. 20, so we will be more specific here. Rather than representing each body 
with just a single polyhedron whose triangular faces form the body's outer surface, as in prior work, we can 
represent each body by multiple polyhedra partially or fully nested inside of one another. Each polyhedron 
has a closed surface consisting of triangular faces, and each triangular face is defined by three vertices. A 
tetrahedron is formed by these three vertices plus the body centroid. This tetrahedron is referred to hereafter 
as a simplex. We require each simplex to have a constant density, but different densities can be assigned 
to each simplex. This allows for density variation over the two angular spherical coordinates within each 
polyhedron, with resolution determined by the size of the faces. The possible nesting of multiple polyhedra 
allows for density variation over the radial spherical coordinate within the body volume as well. It should 
be noted though that to make practical use of the method herein, the coordinates with respect to the body 
centroid of the three non-centroid vertices of each simplex must be known. While this vertex coordinate 
information is known throughout for spacecraft bodies, and can be derived for exterior faces of natural bodies 
from shape data obtained via optical and radar observations or LIDAR surveying, it is difficult to specify 
for interior faces within natural bodies without specific knowledge of internal mass distribution. 

The underlying principle of what follows is that the evaluation of the mutual potential's double volume 
integral over both bodies is equivalent to a global sum of the results of evaluating that double volume 
integral over each possible pairing of simplices, one being drawn from each body. For each such pairing, 
the contribution to the mutual potential is given by a Legendre-series expansion. Successive terms of this 
expansion are linear combinations of other terms which factor into symmetric tensors of increasing rank that 
are independent of relative position and attitude, and other tensors that depend on the relative position and 
attitude. The gradients of the mutual potential according to this formulation make use of efficient tensor 
differentiation rules and the chain rule, resulting in a series expansion for the forces and moments that 
converges for all points exterior to all polyhedra representing the bodies, and hence for all points exterior to 
the bodies. 

It is useful to note here that the number of terms kept in the Legendre-series expansion determines the 
order, in the inverse of the distance between body centroids, of the errors in the forces and the moments. 
When the rigid bodies have small separation distances, more terms are required to maintain the same error 
levels. A possibility for adaptivity during simulation is to adjust the number of terms used in the series 
expansion for the force and moment computations, depending on the body separation distance. Another 
possibility is to refine one or more of the polyhedra representing each body as separation distance decreases. 
However, this must be done carefully, in a manner which involves no change to the relevant conserved 
quantities and no discontinuity in the motion states. 

It is also worth mentioning that one can consider our approach as a special case of a more general 
problem, not treated in this paper, in which the two bodies are no longer rigid but fragmented, i.e. they 
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are so-called "rubble-piles" . The rubble-pile model of asteroids is being supported by an increasing body 
of evidence from recent asteroid exploration missions. 27 A number of researchers have represented such 
rubble-pile asteroids as collections of nonintersecting spheres held together under self-gravity, 28 not to be 
confused with representing a rigid body by filling its volume with spheres fixed with respect to one another. 
One could just as well represent rubble-pile asteroids as collections of ellipsoids or polyhedra of arbitrary size 
and shape held together under self-gravity, rather than smooth spheres, as in Refs. 29 and 30 respectively. 
This allows for filling a body's volume more fully or at different porosities than are possible with spheres 
(even zero initial porosity with polyhedra). With polyhedra, different levels of rigidity within a body can 
then be considered simply by initial combination of smaller polyhedra with coincident or adjacent faces into 
larger polyhedra. The full two rigid body problem that is the subject of this paper can then be viewed from 
a different perspective as the limiting case of that process of combination within two separated bodies. This 
is similar to using a tree-code method, 31 but only using the root- or top-level cell (or group) within each 
body. In other cases employing multiple non-intersecting polyhedra within each body, the gravitational force 
and moment couples between every possible pair of polyhedra in the binary system can still be accurately 
evaluated using the methodology of this section, considered independently from the rest of this paper. 



B. Mutual Gravitational Potential 

We label the two polyhedral rigid bodies B\ and Bi . Consistent with the earlier discussion and definition 
of the polyhedral modelling approach, let body B\ be divided into a set of simplices indexed by a and let 
body &2 be divided into a set of simplices indexed by b . Evaluating the double volume integrals over B\ 
and B>2 is equivalent to the double summation over all a and over all b of the result of evaluating the double 
volume integrals over each simplex combination (a, b). This is shown in the following expression for the 
mutual potential. 18 Note that at this point, we make use of tensor notation and the Einstein convention of 
summation over repeated indices: 
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The scalars T a and TJ, and the tensors Q of increasing rank are all independent of both relative position 
between centroids and relative attitude between the bodies, so they can be computed before any dynamic 
simulation. However the vector w is dependent on relative position, and both the vector w and the matrix 
r are dependent on relative attitude through a matrix v , where 



The Q's arc defined in Rcf. 18, in which they are written out explicitly up to the third rank for illustration 
of their form. The series within the braces in Eq. JQ) is infinite but sufficient accuracy seems to be obtained 
with just the first several terms in square brackets. We denote these bracketed terms as scalars Uo, U\, U2, 
and so on. 



C. Mutual Gravitational Forces and Moments 

The gravitational force existing between the two bodies is given by the partial derivative of the mutual poten- 
tial with respect to the relative position vector X. To obtain this we first derive simple tensor differentiation 
rules: 20 



4 of El 

American Institute of Aeronautics and Astronautics 



These are used in finding each successive dU / dXg term. The first few of these terms are illustrated below: 
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Such terms are used in the overall expression for the force vector that is in turn used in the relative equations 
of motion, namely (laying aside the tensor notation for the moment) 

ax = ~ G ^ Lfl» r -^l ax ax ax -)■ (3) 

The torque or moment between the bodies due to their mutual gravitational interaction is given by an 
expression requiring the partial derivative of the mutual potential with respect to the relative attitude R . 
Evaluating this requires partial differentiation of one attitude rotation matrix, an element of SO (3), with 
respect to another such matrix. For this we use the basic rule 

dRjk _ x4> x k ( 4 ) 



g R 4,e 
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inside of the expression for the tensor D . The details of this tensor and an alternate rule are discussed in 
Ref. 20. It is in turn used in the two tensor differentiation rules: 20 

_ dr^ _ 2 t D « f51 
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These are used in finding each successive dU / dFft term. There is no attitude dependence in Uq, and then 
the next few partials are: 
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Such terms are used in the overall expression 
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Using this within Eq. I|17|) below, we can evaluate the moment due to the mutual gravitational potential. 



III. Full Two Body Dynamics and Continuous Equations of Motion 

In this section, we describe the full two rigid body problem, and we present the continuous relative 
equations of motion (EOM) and the conserved quantities for the full two rigid body dynamics. Continuous 
EOM for the full two rigid body problem are given in Ref. 6, and they are formally derived in the context 
of Lagrangian mechanics in Ref. 25. A more detailed description and derivation of materials in this section, 
including derivation of inertial EOM in addition to relative EOM, can be found in the latter reference. 
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The physical configuration of the full two rigid body problem is described in terms of the position vector 
of each rigid body in an inertial reference frame, an element of R 3 , and the attitude of each rigid body with 
respect to that frame. This attitude is represented by a rotation matrix, which is a 3 x 3 orthogonal matrix 
with determinant +1, and hence is an element of the Lie group SO(3). The configuration space of each rigid 
body is therefore a Lie group referred to as the Euclidean Lie group, SE(3) = R 3 (s)SO(3) (with (s) meaning 
the semi-direct product). The motion of the two rigid bodies depends only on the relative positions and 
the relative attitudes of the bodies. This is a consequence of the property that the mutual gravitational 
potential depends only on these relative variables. Thus, the Lagrangian of the two rigid bodies is invariant 
under a group action of SE(3), and it is natural to reduce the EOM for the full body problem by writing 
them in one of the body fixed frames. Without loss of generality, we present the EOM in the frame fixed to 
B 2 . Then the configuration space of the reduced system is SE(3). 

We start with the relative variables 



X = Ri {XX-X2), 



R — R^Rx, 



(8) 



where X £ R 3 is the relative position of B\ with respect to B 2 expressed in the frame fixed to B2, and 
R G SO (3) is the relative attitude of B\ with respect to B 2 . The corresponding linear and angular velocities 
are 
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where V £ R 3 represents the relative velocity of B\ with respect to B2 expressed in the frame fixed to B2, 
and fl £ R 3 is the angular velocity of B\ expressed in the frame fixed to B 2 . The standard and nonstandard 
moment of inertia matrices (see Rcf. 25 for the distinction) of B\ can also be expressed in the frame fixed 
to B 2 according to Jr = RJ\R T and JdR = RJd x R T ■ Note that Jr and JdR are not constant matrices. 

Assuming that the total linear momentum of the whole system is zero, the kinetic energy of the system 
is expressed in terms of the relative variables as 



KE=±m\\V\\ 2 
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r ^ 1 ^ € R , and S(-) : R 3 1— > so(3) is an isomorphism between R 3 and skew-symmetric matrices, 
defined such that S(x)y — x x y for any x, y £ R 3 . We also have from Eq. 0) the potential energy expressed 
in terms of the relative variables, that is to say PE = U(X, R), recalling that in Eq. r and w are functions 
of X and both w and r are functions of R. The following continuous EOM of the full two rigid body problem 
in relative coordinates can be obtained by taking variations of the reduced Lagrangian, L = KE — PE, as 
follows: 
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where the vector JrQ = RJifli S R 3 is the angular momentum of the first body expressed in the second body 
fixed frame. The moment due to the gravity potential, Af € I 3 , is determined by the following relationship 
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These equations are also given in different notation in Rcf. 20. The total energy, the total linear momentum 
expressed in the inertial frame, and the total angular momentum expressed in the incrtial frame are conserved 
quantities, given by 



7t = m\V\ + m 2 u 2 , 

7tt = x\ x m\v\ + R\JiQ,\ +12 x ra^vi + RaJ2^2- 



E 



KE + PE 



(18) 
(19) 
(20) 



IV. Discrete Equations of Motion: the Lie Group Variational Integrator 

This section presents the LGVI, essentially a statement of discrete relative EOM in Hamiltonian form, as 
opposed to continuous relative EOM in Lagrangian form in the previous section. Again, a detailed description 
and derivation of materials in this section, including development of inertial discrete EOM in addition to 
these relative discrete EOM, can be found in Ref. 25. Here we simply state the results. 

A. Discrete Equations of Motion 

Since the dynamics of the full two rigid bodies arise from Lagrangian or Hamiltonian mechanics, they 
are characterized by symplectic, momentum and energy preserving properties. These geometric features 
determine the qualitative behavior of the dynamics, and they can serve as a basis for further theoretical 
study of the full two rigid body problem. Furthermore, the configuration space has a group structure 
denoted by SE(3). 

However, general numerical integration methods, including the popular Runge-Kutta methods, fail to 
preserve these geometric characteristics. General integration methods arc obtained by approximating con- 
tinuous EOM by directly discretizing them with respect to time. With each integration step, the updates 
involve additive operations, so that the underlying Lie group structure is not necessarily preserved as time 
progresses. This is caused by the fact that the Euclidean Lie group is not closed under addition. 

For example, if we use a Runge-Kutta method for numerical integration of (|13|) . then the rotation matri- 
ces drift from the orthogonal rotation group, SO (3); the quantity R T R drifts from the identity matrix. Then 
the attitudes of the rigid bodies cannot be determined accurately, resulting in significant errors in the gravi- 
tational force and moment computations that depend on the attitude, and consequently errors in the entire 
simulation. It is often proposed to parameterize 1)13(1 by Euler angles or unit quaternions. However, Euler 
angles are not global expressions of the attitude since they have associated singularities. Unit quaternions 
do not exhibit singularities, but are constrained to lie on the unit three-sphere § 3 , and general numerical 
integration methods do not preserve the unit length constraint. Therefore, quaternions lead to the same 
numerical drift problem. Re-normalizing the quaternion vector at each step tends to break the conservation 
properties. Furthermore, unit quaternions double cover SO(3), so that there are inevitable ambiguities in 
expressing the attitude. 

In contrast, the Lie Group Variational Integrator has desirable properties such as symplecticity, mo- 
mentum preservation, and good energy stability for exponentially long time periods. It also preserves the 
Euclidian Lie group structure without the use of local charts, reprojection, or constraints. The LGVI is 
obtained by discretizing Hamilton's principle; the velocity phase space of the continuous Lagrangian is re- 
placed by discrete variables, and a discrete Lagrangian is chosen such that it approximates a segment of the 
action integral. Taking the variation of the resulting action sum, we obtain discrete EOM referred to as a 
variational integrator. Since the discrete variables are updated by Lie group operations, the group structure 
is preserved. Here we present the resulting discrete EOM as follows; the detailed development can be found 



in Rcf. 25. 
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To propagate these equations, we start with a set of initial states, (X , V , R ,il ,Q,2 ), and perform one 
initial evaluation of the mutual potential gradients, obtaining dU /dX and M with Eqs. and 1)170. We 
then find from Eq. 1|21[) . Solving the implicit equations l|22|) and 1|23|) yields the matrix- multiplication 
update matrices F and i 7 ^ for the attitude rotation matrices, and R 1 follows from Eq. (|24|l . After that, 
we use X 1 and R 1 in a new evaluation of the mutual potential gradients. We then compute V 1 , i7 1 , and 
Q 2l from equations f2"5|) . (|2"r)|) and (|2"7j> . respectively. This yields a discrete map (X , V , R , f2 , f2 2o ) i— ► 
(A^ , V^, ijj, fij, f2 2l ), and this process can be repeated for each time step. Note that only one new evaluation 
of the potential gradients is required per time step. The discrete trajectory in reduced variables can be used 
to reconstruct the inertial motion of the bodies. Either concurrently with that propagation or later after 
completion of it, through storing values, we can use the gradient dU/dX, the relative attitude R, and the 
update matrix F 2 with these equations: 

h 2 dU 

^=^+hv 2n + —R n ^-, (28) 

h dU n h dU n 
2^ n djT + 2^ Rn+1 dX~ n 

i? 2 „ +1 =Ri n F 2n . (30) 



+ —Kl^- + —R^^F^, (29) 



In the discrete map defined by the LGVI above, the only implicit parts are Eqs. 1221) and (|23|l . These two 
equations have the same structure, which suggests a specific computational approach. Using the Rodrigues 
formula, we rewrite those equations as an equivalent vector equations, and we solve them numerically using 
Newton's iteration. Numerical simulations show that two or three iterations are sufficient to achieve a 
tolerance ofe = 10~ 15 . 



B. Properties of the Lie Group Variational Integrator 

Since the LGVI is obtained by discrctizing Hamilton's principle, it is symplectic and preserves the structure 
of the configuration space, SE(3), as well as the relevant geometric features of the full two rigid body problem 
dynamics represented by the conserved first integrals of total angular momentum and total energy. The total 
energy oscillates around its initial value with small bounds on a comparatively short timescale, but there is 
no tendency for the mean of the oscillation in the total energy to drift (increase or decrease) from the initial 
value for exponentially long time. In contrast, the total energy behavior with general numerical methods 
such as the Runge-Kutta schemes tends to drift dramatically over exponentially long time. 

The LGVI preserves the group structure. By using the given computational approach, the matrices F n 
and F 2n , representing the change in relative attitude and attitude of B 2 over a time step, are guaranteed to 
be rotation matrices. The group operation of the Lie group SO (3) is matrix multiplication. Hence rotation 
matrices R n and R 2n are updated by the group operation in Eqs. I|24|) and 13U|) . so that they evolve on 
SO(3) automatically without constraints or rcprojection. Therefore, the orthogonal structure of the rotation 
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matrices is preserved, and the attitude of each rigid body is determined accurately and globally without the 
need to use local charts (paramctcrizations) such as Euler angles or quaternions. 

This geometrically exact numerical integration method yields a highly efficient and accurate computa- 
tional algorithm, especially for the full two rigid body problem examined here. In the full two rigid body 
problem there is a large burden in computing the mutual gravitational force and moment for arbitrary bod- 
ies, so the number of force and moment evaluations should be minimized. We have seen that the LGVI 
requires only one such evaluation per integration step, the minimum number of evaluations consistent with 
the presented LGVI having second order accuracy (because it is a self-adjoint method). Within the LGVI, 
two implicit equations must be solved at each time step to determine the matrix-multiplication updates for R 
and i?2- However the LGVI is only weakly implicit in the sense that the iteration for each implicit equation 
is independent of the much more costly gravitational force and moment computation. The computational 
load to solve each implicit equation is comparatively negligible; only two or three iterations are required. 
Altogether, the entire method could be considered "almost explicit" . 

The LGVI is a fixed step size integrator, but all of the properties above are independent of the step size. 
Consequently, we can achieve the same level of accuracy while choosing a larger step size as compared to 
other numerical integrators of the same order. All of these features are revealed in the simulation results 
below. 

V. Numerical Examples 

A. Implementation Details 

Here we describe our specific implementation of the combined computational methods outlined in this paper. 
Our codes consist of pre-processing scripts and post-processing scripts written in the MATLAB scripting 
language, and actual executables written in the C language. 

Body models for a natural body or for a spacecraft are originally provided to us already in the form of 
ordered vertex and face lists. Each row of the vertex list contains three numbers for the coordinates of that 
vertex's position with respect to a specified reference frame, preferably with origin at the center of mass of 
the body and axes aligned with the body's principle axis. Each row of the face list contains the three row 
numbers for the rows of the vertex list corresponding to the three vertices that form the corners of that face. 
The row numbers reference the vertex information in order, moving counterclockwise around the vertices as 
viewed from outside of the body (faces are oriented with outward-pointing normal vectors according to the 
right hand rule). If enough information about density distribution for the body is known, so that a density 
number is also assigned to the simplex associated with each face, those numbers are present in an additional 
column in the ordered face list. Otherwise, all simplices may be given the same density value. In this paper 
we do not address the process of generating such formatted body model files from actual observation data 
for asteroids or from actual CAD software models for spacecraft. The face and vertex lists can be manually 
generated from scratch for simple arbitrary polyhedron shapes, as for the octahedra bodies used in the next 
section. In addition to the body models themselves, the initial conditions are needed, including the initial 
attitude of each body with respect to some common reference frame, the initial spin axis orientation and spin 
rate of each body in that frame, and the initial mutual orbital elements or equivalent relative translational 
motion parameters with respect to that frame. 

The pre-processing scripts handle a number of preliminary computations. For each body, the centroid 
(center of mass) and the principal axis directions are found, and if not at the origin and along the elementary 
unit vectors, respectively, of the frame in which the vertex coordinates are given, the coordinates of every 
vertex are shifted and rotated so that this becomes the case. The moment of inertia contribution of each 
simplex in the new frame is found, and these are summed to obtain the total standard moment of inertia 
matrix. Other information such as the surface area, volume, and mean equivalent radius of each body is also 
found and reported. All vertex coordinates and other quantities can then optionally be nondimensionalized 
using user-entered length, mass, and time factors, for better numerical conditioning. Finally, five files are 
produced. The first two of these files, one for each body, contains a re-ordering of the elements of the position 
vectors, with respect to the centroid, for the vertices of the simplices. The other three files, respectively, 
contain the initial state vector (V , mV , J2^2 j Jro^oi > ^2 ); other system physical data (densities, body 
volumes, J±, J2, mi, m 2 , m, nondimensionalized G, and nondimensionalization factors), and the integration 
parameters (starting and stopping times and truncation error tolerance, if needed). 

A first executable makes use of the Runge-Kutta-Fehlberg 7(8) integration method, hereafter referred to 
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as RKF7(8) or just RKF, to propagate the system, after a starting calculation of the time-invariant Q tensors. 
It is noted that this one-time only calculation of the Q's is analogous to finding the successively higher-order 
mass moments of a "normalized" simplex with vertices at (0,0,0), (1,0,0), (0,1,0), and (0,0,1). It is also 
noted that for this high-order scheme, the EOM are evaluated thirteen times within each integration step, 
and an evaluation of the mutual gravitational potential force and moment is required each of those times. 
After this the state update is performed and the new state vector is written to an output file only if the 
truncation error is within the tolerance specified. Step size adjustment is performed with every step. With 
each state update the mutual potential itself, force, and moment, while not required for the propagation of 
the dynamics, are also evaluated again using the new state vector and written to another output file. This 
allows for checking the total energy conservation and the linear and angular accelerations. Hence fourteen 
force and moment evaluations are needed per integration step. 

Another executable uses the LGVI. At the start of this, the initial state vector from the input file, 
(X ,mV , J2^2 0J Jro^o'-R-o'-R-^o): ^ s converted to the vector for the discrete mapping, (X , V , R , fl , O 2o ). 
Again, rather than fourteen force and moment evaluations, this LGVI involves only one such evaluation per 
integration step, plus the quick solution of the implicit equations. Step size is fixed, but this does not present 
a significant problem for most scenarios observed in binary asteroids having mutual orbits with relatively 
low-eccentricity. 

With the aim of completing simulations much faster than in a single-processor environment, a parallelized 
version of each executable above was also written in C with the addition of MPI. In using the methods 
presented in this paper, most of the computation time is associated with the evaluation of the potential 
gradients, and that involves performing the same operations for all of the different simplex combinations, 
followed by a global sum. This is well-suited for parallelization. The parallelized version of each executable is 
flexible in that any number of nodes or processors can be specified by the user. Then the process assigns to 
each of the other processors the task of calculating the portion of the potential gradients double summations 
of Eqs. and © that arises from pairing a single simplex in £>2 with all simpliccs in B\ in succession. If 
the number of other processors specified by the user or found available on the cluster is less than the number 
of simplices in B2, this is done in rounds until the portion of the problem matching with every a is obtained. 
The parallelized version of each executable has been used on Myrinet clusters at the Center for Advanced 
Computing (CAC) at the University of Michigan and at the Supercomputing and Visualization Center at 
NASA's Jet Propulsion Laboratory. Though compiler and user environment differences produced markedly 
different capabilities in each cluster environment, eventually a further two orders of magnitude reduction in 
computation time over otherwise identical single-processor runs was achieved with both the VI and RKF7(8) 
schemes. It should be noted however, that for the simulations with results presented in the following section, 
the parallel computing capability was not needed for such small (in number of faces) body models, and was 
not utilized. 

The MATLAB script post-processing of the output files generates all desired plots of various dynamic 
quantities, with optional capability for animation generation. All pre- and post-processing steps and scripts 
are identical regardless of whether the RKF7(8) or LGVI executable is used, and regardless of whether the 
single-processor or parallel version of each is used. 

B. Simulation Results 

Simulation results for two octahedral rigid bodies with eight faces and eight simplices each are given for a 
variety of scenarios. Octahedra are used rather than more complex shapes because they are the simplest 
polyhedral shapes that manifest the coupled dynamics behavior desired in all of the scenarios. For greater 
simplicity, the octahedra are made symmetric about all axes, although they are of different sizes. The extents 
data defining the locations of the corners of each octahedron are given in Tabled as are various physical 
parameters of each octahedron including mass and moment of inertia properties. We present simulation 
results for four scenarios, and the initial condition for each scenario is given in Table El 

Scenario 1 The first scenario presented here is that of short-duration simulation of the two octahedra 
starting from initial conditions matching with a medium eccentricity elliptical mutual orbit. Both the 
RKF7(8) and LGVI integrators are used, with the intent of making a direct comparison between the trajec- 
tories of the configuration variables that result from using each integrator over a short simulation duration. 

Figure H shows the difference between the output of the RKF7(8) and that of the LGVI in components 
of reconstructed inertial position, inertial velocity, and body-frame angular velocity vectors for B2, plus the 
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Table 1. Properties of octahedral body models used in simulations 



Property 


£ 2 


B 


| 


Surface area (m 2 ) 




8.839 




2.002 


Volume (m ) 




1.800 




0.1561 


1 J\ 1 LJ.J. V ■ X ClA-li. LAO I ±_L± f 


0.7546 




0.3340 


Mass (kg) 




4500 




390.3 


Density (kg/m ) 




2500 




2500 


t* (kg-m 2 ) 


1377.0 




9.24 


iy.y (kg-m 2 ) 




814.5 




42.99 


7 2Z (kg-m 2 ) 


1462.5 




44.32 


Extents (m) 


min 


max 


min 


max 


body frame X 


-1.0 


1.0 


-1.0 


1.0 


body frame Y 


-1.5 


1.5 


-l/exp(l) 


l/exp(l) 


body frame Z 


-0.9 


0.9 


-1/tt 


1/tt 



Table 2. Initial Conditions 



Scenario 


Attitude* (deg) 


Body spin** (rad/s) 


Orbital elements (m,deg) OR state vector (m,m/s) 


1 


(100, 9.8, 175) 
(160, -5, 165) 


(0, 0, 5.0 x 10" 5 ) 
(0, 0, 9.2 x 10~ 5 ) 


(o,e ) i,n,w,i/)=(4m l 0.3, 5°, 15°, 60°, 10°) 


2 


(180, 0, 30) 
(270, 0, 30) 


(0, 0, 0.566) 
(0, 0, -0.566) 


X = [0,6,0] T , Vq = [-0.006, 0,0] T 


3 


(-22.6, 5, 180) 
(50.3, 5, -180) 


(0, 0, 1.63 x 10" 4 ) 
(0, 0, 1.55 x 10~ 4 ) 


(o,e,i,n,a;,i/)=(52.9m, 0.942, 5°, 0°, 88.2°, -107.1°) 


4 


(-75, 30, 180) 
(-75, 30, 180) 


(0.007, 0.007, 0.05) 
(-0.003, 0.002, 0.004) 


X = [-0.5, 1.8, 1.1] T , V = [-0.3, -0.1, 0] T 



* 3-1-3 Eulcr sequence for Bi (first line) and B2 (second line). 

* * Components of angular velocity of each body expressed in its own body-fixed frame for B\ (first line) and B2 (second line) . 



difference in body attitude parameters for $ 2 . The corresponding output difference plots for body B\ look 
very similar. The differences in vector components of Figure 1(a) are normalized by the system's semi- major 
axis (a = 4.0 m). The differences in vector components of Figure 
circular velocity (\J fi/a = 2.856 x 10~ 4 m/s), and those of Figure 

10~ 5 



1(b) 



1(c) 



are normalized by the equivalent 
are normalized by the equivalent 



meanmotion (-\//V a3 = 7.141 x 10~ 5 radians/s). To obtain the results compared here, the total number of 
mutual potential derivatives evaluations and actual running time using the RKF7(8) routine were 70014 and 
494 seconds, respectively, while the number of such evaluations and actual running time using the LGVI 
were 70001 and 539 seconds. Therefore the computational effort and resources used were roughly the same 
in each case. All of these results show that the LGVI can be trusted to produce almost exactly the same 
trajectory as a standard RKF7(8) integration routine over short time scales. As the simulation duration 
increases the trajectories from the two integrators begin to diverge. The behavior of integrals of motion and 
appropriate error metrics must then be used to discern which trajectory is to be taken as the "truth" . 

Scenario 2 Another scenario is that of propagation from an initial condition with the bodies aligned but 
possessing relatively large magnitude centroid velocity vectors that are antiparallel and perpendicular to the 
initial line between centroids. This scenario is simulated both with the LGVI at different step sizes and 
with the variable step size RKF7(8) at different error tolerances. This allows for a comparison between the 
integrators of their performance, in terms of the total energy and total angular momentum integrals and 
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4000 0000 
Time (sec) 



4000 6000 
Time (sec) 



(a) Inertial position for B'z 



(b) Inertial velocity for B2 




10000 




2000 4000 6000 8000 

Time (sec) 



10000 



(c) Angular velocity for B2 



(d) 3-1-3 Euler angles for B2 



Figure 1. (Scenario 1) Difference between RKF7(8) and LGVI output. 



the attitude error metric growth vs. computational burden. The results in Table |3] illustrate the general 
superiority of the LGVI approach over Runge-Kutta-type approaches. 

Here we see that for any pair of simulations, one using the LGVI and the other using the RKF7(8) scheme, 
for which the total energy metric performs about the same, the computation time needed to complete 
the simulation using the LGVI is a fraction of that needed using the RKF7(8). Simultaneous with this 
improvement in run time, the total angular momentum and attitude error metrics still perform better in 
the LGVI run than in the RKF7(8) run by multiple orders of magnitude. Going in the other direction, 
as the step size for the LGVI is reduced so that the computational burden using it begins to approach 
that for any chosen run using the RKF7(8), all error metrics remain at the same level as or else orders of 
magnitude smaller than those for the chosen RKF7(8) run. For the LGVI, the round-off error accumulates 
when multiplying rotation matrices at 1|24|) . The rotation matrix error of the LGVI is caused only by the 
floating-point arithmetic operation, and it is increased as the number of integration steps is increased. A 
similar trend is observed in the total angular momentum error for the LGVI, because determination of the 
total angular momentum in the inertial frame from the states written to the output file makes use of the 
rotation matrices. 



Scenario 3 The next scenario illustrates the ability of our methods to capture the interesting effects of 
coupling in a mutual orbit configuration that the Kcplerian two-body approximation incorrectly predicts as 
being perpetual. Simulation with the LGVI yields the trajectory illustrated in Figure|2(a)| which transitions 
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Table 3. (Scenario 2) Performance comparison between RKF7(8) and LGVI 



Method 


h* 


N* 


l w 


e« 


E[|ATE|]tt 


EOIAtpHI] 


+ 


E[ 


I-R T R 




RKF7(8) 


0.236 


2368912 


23439 


10" 12 


3.901 x 10" 


12 


1.493 x 10 


—9 


1.151 x 10" 


1 


RKF7(8) 


0.421 


1331414 


9102 


lO" 10 


1.274 x 10" 


10 


2.630 x 10 


- 7 


1.985 x 10" 


5 


RKF7(8) 


0.749 


747376 


5252 


10~ 8 


2.284 x 10" 


-8 


4.620 x 10 


— 5 


3.173 x 10- 


3 


T PUT 

LKj VI 


U.Uloy 




1 QC1 1 




1.698 x 10- 


11 


0.10/ X 1U 


-10 


2.525 x 10- 11 


LGVI 


0.04 


1000000 


9920 




1.928 x 10- 


11 


1.189 x 10" 


-10 


2.120 x 10~ n 


LGVI 


0.08 


500000 


5127 




9.879 x 10- 


11 


4.139 x 10" 


-11 


2.004 x 10~ 12 


LGVI 


0.4 


100000 


983 




2.234 x 10" 


-9 


6.266 x 10" 


-12 


3.386 x 10~ 14 


LGVI 


0.8 


50000 


431 




9.326 x 10" 


-9 


1.279 x 10" 


-11 


6.352 x 10~ 14 


LGVI 


1.0 


40000 


335 




1.512 x 10" 


-8 


3.991 x 10" 


-12 


4.786 x 10~ 14 



* h is integration step size, in seconds, fixed for LGVI but averaged over the run's duration for RKF7(8) 

* N u is the total number of calculations of the mutual potential derivatives made within the run 
°tw is the "wall-clock" time to complete each simulation run, in seconds 

" e is the error tolerance for the variable step size in RKF7(8) 

t TE and ttt are total energy and the total angular momentum, respectively, while A refers to deviation from the 

initial value over simulation 
t E[-] denotes mean 



from a highly elliptical orbit to a hyperbolic escape path. This is shown by the plots in Figure |2(b)| of the 
semi-major axis and eccentricity change during the close encounter, which occurs roughly midway through 
the run duration of 60,000 seconds. The initial conditions and body configurations arc symmetric about the 
initial orbital plane, and as such the motion of the centroids should be restricted to the initial orbital plane. 
This is observed numerically, as the body centroids remain within 8.6 x 10 -14 meters of the initial orbital 
plane throughout the simulation. 




(a) Trajectory of binary octahedra system components (b) Eccentricity and semi major axis 

Figure 2. (Scenario 3) Disruption of the binary octahedra system. 

Scenario 4 Finally, we examine a very long duration simulation starting from initial conditions that are 
stable in the sense that orbital trajectories are confined to separated and bounded regions but are highly 
unstable in the sense that body attitudes vary greatly and irregularly. The trajectory of binary octahedra 
system is shown in Figure [3] Note that this motion cannot be observed with the point mass assumption 
in the classical two body problem. This scenario is simulated both with LGVI at different step sizes and 
with RKF7(8) at different error tolerances, for 5 x 10 6 (sec) maneuver time. Table 0] illustrates the general 
superiority of the LGVI approach, similar to Scenario 2. In addition, we see that the simulation results of 
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Y(m) " 10 "5 X (m) 

Figure 3. (Scenario 4) Trajectory of binary octahedra system 



Table 4. (Scenario 4) Performance comparison between RKF7(8) and LGVI 



Method 


h 


N u 


t w 


e 




E[ ATE ] 




E[||A7rr||] 




E[ 


I-R T R 




RKF7(8) 


8.81 


7373106 


48850 


10" 


16 


6.198 x 10 


-8 


3.830 x 10" 


-6 


6.089 x 10- 5 


RKF7(8) 


20.66 


3145532 


20873 


10" 


13 


8.731 x 10~ 


-5 


7.427 x 10' 


-3 


1.015 x 1Q- 1 


RKF7(8) 


48.74 


1333475 


9699 


io- 


10 


1.959 x 10" 


-1 


5.246 x 10 





1.834 x 10 1 




LGVI 


0.7 


7142858 


46963 






6.947 x 10- 


11 


2.190 x 10- 


12 


6.645 x 10- 13 


LGVI 


2 


2500000 


17320 






5.687 x 10- 


10 


5.077 x 10- 


13 


3.663 x 10- 13 


LGVI 


5 


1000000 


6897 






3.601 x 10~ 


-9 


8.447 x 10- 


13 


1.642 x 10- 13 


LGVI 


10 


500000 


3750 






1.517 x 10" 


-8 


2.567 x 10- 


13 


3.521 x 10" 13 



See Table 131 for notations. 



RKF7(8) with larger error tolerances, 10 -13 and 10 -10 , are completely unreliable since the rotation matrix 
error is increased to the unacceptable levels of 10 _1 and 10 1 . 

The quantitative comparison is summarized in Figure where the total computation time and the stan- 
dard deviation of the total energy are shown over the number of the mutual potential derivatives calculations. 
In Figure |4(a)| we see that the total computation time is directly proportional to the number of the mu- 
tual potential derivative calculations regardless of the integration scheme used. This justifies the statement 
that the LGVI is an "almost explicit" computational method; the computational load to solve the implicit 
equations is comparatively negligible. Considering that the mutual potential derivative computations are 
the major computational burden for the numerical simulation of full two rigid body dynamics, an efficient 
integration scheme with fewer mutual potential derivative calculations should exhibit smaller error measures. 
In Figure [4(b)| the lower left portion represents higher computational efficiency, and the upper right portion 
represents lower efficiency. We see that LGVI is more efficient than RKF7(8). Figure compares the total 
energy deviation history and the relative rotation matrix error history. For LGVI, the total energy varies 
within the level of 10 -11 . The repeated peaks of the total energy deviation correspond to the perigees of the 
orbit shown in Figure |3 but but there is no tendency for the mean of the variation to drift for the entire 
simulation time. The rotation matrix error is slightly increasing, since the round off error in multiplying 
orthogornal matrices at (|24(l is accumulated. But the error is below 10 -12 over the entire simulation's span. 
However, for RKF7(8), the variation of the total energy is linearly increasing over time, and the rotation 
matrix error is increasing to the level of 10" 4 . We have seen that LGVI exhibits good total energy behavior 
for exponentially long time periods, and it also preserves the group structure well. The accuracy of RKF7(8) 
is vulnerable to a long time simulation of the full two rigid body dynamics. 
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■ RK7(8) 

• LGVI 



012345678 
Total number of potential derivative calculation x ^ 



V 

d 

& 

1 10"' 
o 



RK7(8) 
LGVI 



012345678 
Total number of potential derivative calculation x jq 6 



(a) Running time v.s. number mutual 
potential derivatives calculations 



(b) Mean deviation of TE v.s. number of 
mutual potential derivatives calculations 



Figure 4. (Scenario 4) Quantitative comparisons of running time and standard deviation of total energy over 
number of mutual potential derivatives calculations. 



.x 10 



.x 10' 



1.5 



m i 



) 

x 10~ 12 


1 2 


3 4 

x ID 6 















2 3 
time (sec) 




x 10 



2 3 
time (sec) 



x 10 



(a) LGVI 



(b) RKF7(8) 



Figure 5. (Scenario 4) Qualitative comparisons of behavior of deviation in total energy and rotation matrix 
error over time with similar computational load for LGVI and RKF7(8). 



VI. Conclusions 



This paper presents a complementary combination of tools or algorithms that achieves superior perfor- 
mance in the computationally demanding simulation of the full two rigid bodies problem: a simple method 
for efficient calculation of the mutual gravitational potential and its derivatives given versatile polyhedral 
models of the rigid bodies and a Lie group variational integrator consisting of discrete relative equations 
of motion that preserves the geometric features and structure of the configuration space. The use of these 
easily implemented methods together allows for simulation of the full two rigid bodies, capturing their fully 
coupled dynamics, that is both accurate and efficient. The results above show the accuracy maintained for 
energy, momenta, and attitude geometry constraints, even for long simulation run times. The presented 
mutual potential computations for polyhedral approximations are themselves efficient, but for bodies of rel- 
evant (i.e. large) model size, they still represent the bulk of the computational burden during simulation, 
regardless of the integration scheme used. Hence use of the Lie Group Variational Integrator, which requires 
minimal mutual potential computations, achieves an even greater combined efficiency. This yields maximal 
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simulation durations for fixed computational resources. 
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